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Abstract 

A fractional step method for the solution of steady 
and unsteady incompressible Navier- Stokes equations 
is outlined. The method is based on a finite-volume 
formulation and uses the pressure in the cell center 
and the mass fluxes across the faces of each cell as 
dependent variables. Implicit treatment of convective 
and viscous terms in the momentum equations enables 
the numerical stability restrictions to be relaxed. The 
linearization error in the implicit solution of momen- 
tum equations is reduced by using three subiterations 
in order to achieve second order temporal accuracy for 
time-accurate calculations. In spatial discretizations 
of the momentum equations, a high-order (3 rd and 
5 th ) flux-difference splitting for the convective terms 
and a second-order central difference for the viscous 
terms are used. The resulting algebraic equations are 
solved with a line-relaxation scheme which allows the 
use of large time step. A four color ZEBRA scheme 
is employed after the line-relaxation procedure in the 
solution of the Poisson equation for pressure. This 
procedure is applied to a Couette flow problem us- 
ing a distorted computational grid to show that the 
method minimizes grid effects. Additional benchmark 
cases include the unsteady laminar flow over a cir- 
cular cylinder for Reynolds Numbers of 200, and a 
3-D, steady, turbulent wingtip vortex wake propaga- 
tion study. The solution algorithm does a very good 
job in resolving the vortex core when 5 <h -order up- 
wind differencing and a modified production term in 
the Baldwin-Barth one-equation turbulence model are 
used with adequate grid resolution. 
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1. Introduction 

Numerical solution of the Incompressible Navier- 
Stokes (INS) equations requires special attention in 
order to satisfy the divergence-free constraint on the 
velocity field. One way to avoid the numerical diffi- 
culty is to use an artificial compressibility method. 1-5 
With the artificial compressibility method, the elliptic- 
parabolic type equations are transformed into hyperbolic- 
parabolic type equations. Well-established solution al- 
gorithms developed for compressible flows can be uti- 
lized to solve the resulting equations. In addition, 
an unstaggered grid orientation can easily be incorpo- 
rated in the artificial compressibility approach making 
the method simple and robust. Time-accurate solution 
of INS equations with artificial compressibility method 
requires, in a sense, the solution of a steady-state prob- 
lem in order to advance one physical time step. This 
is possible with subiterating the equations in pseudo- 
time until divergence-free constraint on the velocity 
field is satisfied at each physical time step. When one 
needs to use small physical time step in time-accurate 
calculations in order to capture the details of the flow 
physics, the artificial compressibility approach can be 
costlier. In order to reduce the number of subitera- 
tions, one needs to experiment for the optimum range 
of the artificial compressibility coefficient, 0 . In the- 
ory, the artificial compressibility approach should give 
the same solution at the steady-state no matter what 
value of 0 is used. However, the value of 0 affects 
the convergence rate of the algorithm. Use of implicit 
and Newton-like iteration schemes reduces this prob- 
lem substantially. 5 Further improvement in the con- 
vergence rate of artificial compressibility method for 
unsteady flows is an ongoing research effort. 

An alternative to the artificial compressibilty 
method is a fractional step method, which is espe- 
cially suitable for time-accurate calculations. In the 
fractional step method, the auxiliary velocity field is 
obtained by solving the momentum equations. Then, 
a Poisson equation for pressure is formed by taking 
the divergence of the momentum equations and by 
using a divergence-free velocity field constraint. The 
numerical solution of the Poisson equation for pres- 
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sure with the Neumann-type boundary conditions ex- 
ists only if the compatibility condition is satisfied. In 
three-dimensional curvilinear coordinates, solving the 
resulting algebraic equations from Poisson and mo- 
mentum equations efficiently is one of the important 
features of the fractional step method. Using staggered 
grid orientation might also be neccessary to satisfy the 
continuity equation. 6-13 

The staggered formulation of the fractional-step 
method has been successully applied in Cartesian 
coordinates. 7 , but in a body-fitted curvilinear coor- 
dinates, the formulation is not straight-forward. Shvy 
et al. 8 attempted to develop a method for general ge- 
ometry using contravariant-type velocity components 
in the correction step. Rosenfeld et a/. 9-10 success- 
fully developed a fractional-step solution procedure 
in curvilinear coordinates using the pressure and the 
mass fluxes, divided by constant density, across the 
faces of each cell as dependent variables. The re- 
sulting solver was validated using smoothly varying 
and nearly orthogonal grids. This solver was intented 
to be used for time-dependent flows requiring small 
physical time steps. However, there are CFL num- 
ber restrictions for a wide variety of applications, and 
the method requires carefully generated grids without 
metric discontinuities and with very small streching ra- 
tios. In the present paper, the idea of solving the mass 
fluxes as dependent variables in momentum equations 
is used and the stability of the numerical method is 
greatly improved by treating both convective and vis- 
cous fluxes implicitly. The factorization error is re- 
moved by using a relaxation scheme to solve the al- 
gebraic equations. The time integration scheme and 
operator splitting technique are formulated such that 
overall second order temporal accuracy is mantained 
by reducing linearization error in the implicit proce- 
dure. The validation cases are selected to demonstrate 
that the solution procedure is robust under severe grid 
conditions with metric discontinuity. It is shown that 
large CFL numbers can be used for steady-state as 
well as time-accurate solutions. To demonstrate this 
capability, computed results for a wingtip wake vortex 
flowfield is presented in detail. 

In the following section, the governing equations 
and the present fractional step algorithm for the so- 
lution of the incompressible Navier-Stokes equations 
will be outlined. Next, validation cases in two- and 
three-dimensional space will be presented. 


2. Method of Solution 

For simplicity, fractional-step procedure in sec- 
tion 2.2 is described in Cartesian coordinates. The for- 
mation of conservation of mass and momentum equa- 


tions and their spatial discretizations are outlined in 
section 2.3. Section 2.4 details the turbulence model 
used in tip vortex wake propagation study. 


2.1. Governing Equations 

The unsteady, incompressible Navier-Stokes equa- 
tions are comprised of mass conservation, 


dui 

dxi 


= 0 


and momentum conservation equations. 


( 1 ) 


duj _ dp _ diijUj drij_ _ dp 

dt dxi dxj dxj dxi 


( 2 ) 


In this case, t is the time, the Cartesian co- 
ordinates, Ui the corresponding velocity components, 
and p the pressure. For a Newtonian fluid, the viscous 
stress tensor can be written as 


dui duj 

T ” = *' ( ai7 + ai7 ) 

where u is the effective kinematic viscosity. 


(3) 


2.2. Fractional-Step Procedure 

The time integration scheme is based on oper- 
ator splitting, which can be accomplished in several 
ways by combining the pressure, convective, and vis- 
cous terms in the momentum equations. The frac- 
tional step method, a projection method developed by 
Chorin 11 , is based on the decomposition of vector field 
into a divergence free component and gradient of a 
scalar field. The common application of this method 
is done in two steps. The first step is to solve for an 
auxiliary velocity field using the momentum equations. 
In the second step, the velocity field is corrected by us- 
ing the pressure which can map the auxiliary velocity 
onto a divergence free velocity field. The momentum 
equations are discretized in time using a three-point- 
backward difference formula: 

4(3<-K + <-‘) = -fi + RK) (4) 

where u * denotes the auxiliary velocity field. The 
R(u * ) term in the momentum equations includes the 
convective and viscous terms. It should be noted that 
the time derivatives can be differenced using the Euler 
backward formula for steady-state calculations. The 
velocity field which satisfies the incompressibility con- 
dition has been obtained by using the following cor- 
rection step. 
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— u*) = _Vp' (5) 

where p = p n+1 — p n . At n + 1 time level, the veloc- 
ity field has to satisfy the incompressibility condition 
which is the continuity equation. 

V ■ u n+1 = 0 (6) 

This incompressiblity condition is enforced by using a 
Poisson equation for pressure. 


The mass conservation equation is evaluated over the 
faces of a computational primary cell with volume 
which is shown for Pjki in Figure 1-a. Each term in 
equation (10) approximates the mass flux over the cor- 
responding cell face. If the mass fluxes are chosen as 
unknowns, the continuity equation is satisfied auto- 
matically in generalized coordinate systems. The mass 
fluxes over the £, p , and £ faces of the computational 
cell are 

U f =S f • u 


vV = iv.„‘ (7) 

The Poisson equation for pressure is obtained by tak- 
ing the divergence of equation (5) and using equation 
( 6 )- 

In equation (4), both convective and viscous terms 
are treated implicitly. In order to maintain second or- 
der temporal accuracy, the linearization error in the 
implicit solution of equation (4) needs to be reduced. 
This is achieved by using subiterations. In most cases, 
three subiterations are sufficient to reduce the lineariza- 
tion error. It should be noted that the purpose of 
subiteration procedure here is quite different than in 
the artificial compressibility method. Artificial com- 
pressibility formulation requires the solution of a steady- 
state problem at each physical time step. Therefore, 
the number of subiterations in an artificial compress- 
ibility approach can be an order of magnitude higher 
than number of subiterations for the present formula- 
tion. 


2.3. Discretization 


Since the spatial discretization is based on the finite- 
volume formulation, the governing equations will be 
written in the integral form for the conservation of 
mass 


i 


u - dS = 0 


( 8 ) 


and momentum. 



£ (— uu — pi 4- i/(Vu -f (Vu) T )) • dS 
j>f-dS 


( 9 ) 

where S is the surface area vector, V is the cell volume, 
and (.) T is the transpoze operator. The discretization 
of the mass conservation equation (8) in finite volume 
formulation gives 


(S f • - (S f • u)j_i t t ,+ 

(S’ 7 • U )j,t+A,/ — (S’* • u)j- t _i j+ (10) 

( s< • u )i,fc,J+i “ ( s< = 0 


U v =S V ■ u (11) 

U c =S C • u 



Figure 1-a. Schematic of a primary cell and 
staggered grid orientation. 


Jd. k*l/2. 1-172 



Figure 1-b. Computational cell for U* momen- 
tum equation. 
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The continuity equation with this choice of the de- 
pendent variables takes a form identical to the Carte- 
sian case. Therefore, the mass fluxes are considered as 
the ‘natural’ dependent variables for projection meth- 
ods in curvilinear coordinates. The mass conservation 
equation with new dependent variables in a general- 
ized coordinate system becomes 


u U‘.i- u U‘,- +v h+»- u ]*-b' + 






( 12 ) 


Treating the mass fluxes as dependent variables in 
finite volume formulation is equivalent to using con- 
travariant velocity components, scaled by the inverse 
of the transformation Jacobian, in a finite-difference 
formulation. The mass fluxes in finite-volume dis- 
cretization were used in references 9 — 10, and con- 
travariant velocity components in finite-difference dis- 
cretization were used in reference 12. The choice of 
mass fluxes as dependent variables complicates the dis- 
cretization of the momentum equations. In order to 
replace u by the new dependent variables [/', the cor- 
responding area vectors are dotted with the momen- 
tum equations. Then the integral momentum equa- 
tion is evaluated on different computational cells for 
each unknown U ! . Each cell has the dimensions of 
A£ x A x A£ , but the centers are located at (j + 1 , k ,1) , 
( j , k + /), and ( j , k, l + 5) for , U v , and mo- 

mentum equations, respectively. The computational 
cell with volume Vy + 1 / 2 , fc , r for {/^-momentum equation 
is shown in Figure 1-b. The staggered grid orienta- 
tion eliminates pressure checker-board-like oscillations 
in pressure and provides more compact stencils. The 
derivation of {/^-momentum equation will be outlined 
in this section. The U n -, and 17 -momentum equations 
can be obtained by using cyclic permutation. Spa- 
tial discretization of the momentum conservation law 
equation (9) for a computational cell with volume V 
yields 

= J2 s ‘ f ( 13 ) 

The dot product of equation (13) and results in 


d(Vu) 

dt 


d(vut) 

dt 




(14) 


where the summation is over all the faces of a compu- 
tational cell. It is to be noted that 


u = S e ^ + S j; C/' ) -i-S c t/ c (15) 


The invariance of the velocity vector requires: 

S' • S m = 6i m (17) 

where is the Kronecker delta, and S' is the in- 
verse base to S m . A uniform velocity field can be nu- 
merically preserved if the covariant base vector S m is 
computed at each point from the relation 

S m = (S m+1 x S m+2 )/(S™ • (S m+1 x S m+1 )) (18) 

which satisfies (17) identically. The variable m is the 
cyclic permutation of (^, »7, C)- 

In constructing momentum equations, the product S' • T 
should be computed for each face of each momentum 
equation (see equation 14). For example, the £ face- 
center for the momentum cell is located at (j, k, /). 
The flux over this face is computed from 


(S'T^.fc.i = (-C/ft/'s'-Sfp+S«^(Vu+(Vu) T )) jii; , 

( 19) 

The conservative form of the velocity vector gradient 
is 


Vu = 


_ is uds 

v 


( 20 ) 


Applying equation (20) for the computation of Vuj^ j 
yields 


Vu = — (S ? , , , ,u,- . 1 ut — j . ,u,_i l 1 

y v j + ±,k,l J + J s>*>‘ 

+ s J, t+ i ,“,-,*+4.1 - S It.i'Mik-iJ ( 21 ) 


The t] and £ face centers are located at (j 4- 1 / 2 , k — 
1/2,1) and (j+ 1/2, k, l- 1/2), respectively. The fluxes 
over these faces are computed in a similar way. The 
convective and viscous fluxes in equation (19) may be 
approximated in various ways. In the present work, 
the viscous fluxes are computed by simple averaging 
which results in second order central differencing. The 
convective flux terms, U^U'S 1 in equation (19), are 
computed using higher-order upwind-biased stencil. 
Flux-difference splitting is used here to structure the 
differencing stencil based on the sign of the eigenvalues 
of the convective flux Jacobian. The numerical flux 
fj+1/2 f° r the convective term in equation (19) is given 
by 


and 


U l - S' • u = S' • S m U m 


(16) 
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fj+ 1/2 = 2 [/( u i+0 + /( u ;) ~ 1/2] (22) 


where the <f>j + 1 / 2 is a dissipation term. For 4>j + 1/2 
= 0 this represents a second-order central difference 
scheme. A first-order upwind scheme is given by, 

^+i / 2 = [A/ ; + +1/2 -A^- +1/2 ] (23) 

and a third-order upwind flux is defined by 

^i+1/2 = ~~[A //_ l/2 - A/+ 1/2 + 

A/j+ 1/2 — A/j +3 / 2 ] 

A fifth-order-accurate, upwind-biased stencil which re- 
quires only seven points was derived by Rai 14 as follows 

t j+ il2 = -±[-2Af+_ 3/3 + HAtf_ 1/2 

-6A/ i f f i/2 - ZAff +3/2 ^5) 

+2A/ j+5 ^ 2 — HA fj +3 / 2 
+6A/ j - +1 ^ 2 + 3A fj_ 1 / 2 ] 

where A /* is the flux difference across positive or neg- 
ative traveling waves. The flux difference is computed 
as 

A f?+i/2 = ^(^Auj+j^ (26) 

where the A operator is given by 

Au;+i /2 = -Uj+1 - Uj (27) 

The plus (minus) Jacobian is computed from 

a± = \( a ± l a l) ( 28 ) 

The Roe properties 15 which Me necessary for a conser- 
vative scheme, are satisfied if the following averaging 
procedure is employed 

“= |(«j+i + “i) (29) 

An implicit, delta-law form approximation to the mo- 
mentum equations after linearization in time results in 
a hepta-diagonal scalar matrix equation written as 

b£qj-i,k,l + af>9j,k,i + c&qj+i,k,i + d6qj t k-u 
+ef>qj,k+i,i + f6qj,k,i-i + <jf>qj,k,i + 1 = RH.S. 

(3°) 

where 6q = U n+1 — U n and a,b,c,d,e, /, and g are 
diagonals. The Gauss-Seidel line relaxation scheme, 
which was successfully employed by MacCormack 16 , 
is used to solve the matrix equations. In Eq. (30), 
the right-hand-side term is computed and stored for 


the entire domain. The line relaxation procedure is 
composed of three stages, each stage involving a scalar 
tridiagonal inversion in one direction. In the first stage, 
6q is solved line-by-line in one direction. Before the 
tridiagonal equation is solved, off-tridiagonal terms are 
multiplied by the current value of 6q and are shifted 
over to the right-hand-side of the equation. The same 
procedure is repeated in the second and third stage 
by inverting the tridiagonal matrix in one direction, 
and treating the off-diagonal terms for the other two 
directions in Gauss-Seidel fashion. One forward and 
one backward sweep in each computational direction is 
sufficient for most problems, but the number of sweeps 
can be increased. 

A pressure Poisson equation for volume Vj.fc.i is 
formed by taking the divergence of the pressure gradi- 
ent terms in momentum equations (see equation (19) 
and (5)). The resulting algebraic matrix equation in 
generalized curvilinear coordinates, contains nineteen 
diagonals. After two-sweeps of Gauss-Seidel line- 
relaxation procedure in each direction, a four color 
ZEBRA scheme, used in reference 9, is also utilized 
in the present work for the efficient solution of the 
Poisson equation. 

2.4. Turbulence Model 

The turbulent flow calculations use the one- 
equation turbulence model developed by Baldwin and 
Barth. 17 In this model, the transport equation for the 
turbulent Reynolds number is derived from a simpli- 

r ir r ii. . j if. j_i A : 'T'V . « 

lieu lurili ui tut; stdiiueuu k — t muuci ctjUdtiuua. xuc 

model is relatively easy to implement because there 
is no need to define an algebraic length scale. The 
transport equation is solved by using the same Gauss- 
Seidel type line-relaxation scheme as the mean-flow 
equations. The details of this model can be found 
in reference 17. The wake- vortex calculations in sec- 
tion (3.3) use various approximations of the produc- 
tion term in the Baldwin-Barth one-equation model. 
The production term P for vRt is given by 

P ~ cvR t X (31) 

where c is a constant, v is the kinematic viscosity, Rt 
is the turbulent Reynolds number, and X is a scalar 
parameter which needs to be determined. Originally, 
the Baldwin-Barth one equation turbulence model was 
developed by using the magnitude of vorticity as X. 

X = |w| (32) 

Spalart 18 , and Dacles-Mariani 27 suggest to combine 
the magnitude of vorticity |u>| and the strain rate |s| = 
(2 SijSij) 1 ! 2 as follows 
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(33) 


A" = M + xfactor{min{ 0, |s| — |w|)) 

In reference 27, xfactor is set to 2.0 for tip vortex 
calcualtions. Another option for X is to use mimimum 
value of the magnitude of the vorticity and the norm 
of the strain rate tensor 

* = + < 34 > 

The computed results obtained using these production 
terms are compared in section 3.3. 


3. Computed Results 

Computed results for two laminar cases, and one 
turbulent case are presented in this section, Couette 
flow on a distorted grid, unsteady flow over a circu- 
lar cyclinder, and a three-dimensional wingtip vortex 
wake flow are discussed. The wake-vortex calculations 
use the one-equation turbulence model developed by 
Baldwin and Barth. 17 The computer time requirement 
is 50.0 x 10 -6 sec/grid point /iteration on Cray-C90 
computer. Presently, the code has not been fully vec- 
torized, so this figure will be reduced in the future. 

3.1. Laminar Couette Flow 

The first validation case is a grid quality study for 
laminar Couette flow, originally studied in reference 
19. The computational grid with 63 x 63 mesh points 
is shown in Figure 2. The grid is intentionally gener- 
ated in a saw-tooth shape to introduce metric disconti- 
nuity and non-orthogonality. Even with this “not-so- 
nice” mesh point distribution, the scheme should be 
able to resolve linear u- velocity profile which is the 
exact solution for the laminar Couette flow. The flow 
is started with freestream velocity everywhere except 
the stationary wall. The stationary wall has no-slip, 
and upper wall has a moving wall boundary condition. 
Inflow and outflow boundaries are periodic. CFL num- 
ber of 100 is used for this computation, where the CFL 
number is defined as 

CFL = max(\U*\ + \U n \ + |t/ c |) * dt/V (35) 

Figure 3 shows axial ( U ) velocity contours at steady- 
state. The velocity contours show very small kinks 
where metric discontinuties are present in the mesh. 



Figure 2 : Computational grid with 63x63 mesh 
points for Couette flow. 




Figure 3 : U-velocity contours for Couette flow. 


The U— velocity profile at x/L = 0.5 station 
compares very well with the exact solution of the Cou- 
ette flow, as shown in Figure 4. This case shows that 
the current approach introduces minimal grid effects 
where a sudden change in slope of grid lines occurs. 
The small grid quality errors that do arise have an 
insignificant effect on the solution shown in Figure 4. 
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Figure 4 : U-velocity profile at for Couette flow. 



Figure 5 : Couette flow convergence history for 
present fractional-step formulation 

The convergence history of the present method for this 
case is plotted in Figure 5. The solid line shows the 
maximum residual of the momentum equations, and 
the dashed line represents the maximum divergence 
of the velocity. The maximum divergence of velocity 
curve flattens about 10~ 5 value because the iteration 
procedure in the solution of the Poisson equation is 
terminated after achieving e = 10 -5 accuracy. With 
further iterations in the Poisson equation solution, the 
error in the divergence of velocity can be lowered to 
machine zero. Since this is a steady-state case, the 
Poisson equation has not been iterated beyond errors 
below 10~ 5 . It is observed that the convergence behav- 
ior of the existing method is good for this steady-state 
case with large CFL number. In addition, the diver- 
gence of velocity error is guaranteed below an e value 
at each time step in fractional step method. 


3.2. Vortex Shedding From a Circular Cylinder 

The unsteady laminar vortex shedding from a 
circular cylinder at a Reynolds number of 200 was 


studied as a second test problem. The Reynolds num- 
ber is based on the cylinder diameter and the free 
stream velocity. A 129x129 point O-type grid is used 
for this calculation, and the far field boundary is ex- 
tended 15 diameters from the cylinder. The time step 
At of 0.025 w as used and an c value of 10 -5 was set 
as a maximum error tolerance in the solution of Pois- 
son equation. The flow was started impulsively from a 
free-stream condition and run until a periodic vortex 
shedding in the wake occured. The time evaluation 
of the lift and drag coefficients is plotted in Figure 
6. The asymmetric wake started to develop within a 
nondimensional time of 40, and a nearly periodic solu- 
tion is obtained by a nondimensional time of 100. The 
values of the aerodynamic coefficients for the periodic 
cycle are: Cl = 0.0 ± 0.67 and Cd = 1.26 ± 0.04. The 
Strouhal number can be calculated from the frequency 
of the oscillations for the lift coefficient and is found to 
be 0.184. Lift and drag coefficients and Strouhal num- 
ber from the current computations are compared with 
numerical results by Rogers 4 , Rosenfeld 9 , Lecointe and 
Piquet 20 , Martinez 21 , Lin 22 , Thoman and Szewczyk 23 , 
and with experimental data by Experimental values by 
Wille 24 , Kovasznay 25 , and Roshko 26 in Table 1. The 
values of Strouhal number for computed results and 
experimental data are in the range of 0.160 and 0.227. 
Numerical results obtained from the current approach 
are closest to the artificial compressibilty results of 
Rogers 4 (5th order) even though the two methods are 
formulated quite differently. The present results also 
agree well with experimental data 24-26 . 



Figure 6 : Lift and drag coefficients versus time 
for flow over a circular cylinder at a Reynolds 
number of 200. 

In Figure 7, particle traces and vorticity magni- 
tude contours are plotted at various time during one 
period, which occurs over a nondimensional time of 
5.43. 
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Cl 

Cd 

St 

Present 

±0.67 

1.26 ±0.04 

0.184 

Rogers (num) 
3rd order 
5th order 

+0.75 

±0.65 

1.29 ±0.05 
1.23 ±0.05 

0.160 

0.185 

Rosenfeld (num) 

+0.70 

1.40 ±0.04 

0.201 

Lecointe & Piquet (num) 
2nd order 
4th order 

±0.70 

±0.50 

1.35 ±0.04 
1.58 + 0.0035 

0.227 

0.194 

Martinez (num) 


1.27 ± 0.0035 


Thoman & Szewczyk (num) 

1.17 + 0.005 


Witle (exp) 
Kovasznay (exp) 


1.30 

0.19 

Roshko (exp) 



0.19 


Table 1 : Lift and drag coefficients for flow over 
a circular cylinder at a Reynolds number of 200. 



Figure 7 : Particle traces and vorticity mag- 
nitude contours for flow over circular cylinder 
at a Reynolds number of 200 at various nondi- 
mensional times. 

The first plot shows the extension of the top vortex, 
resulting in maximum values of lift and drag. The next 
plot corresponds approximately to mimimum in drag 
with zero lift. The third and fourth plots are mirror 
images of the first two plots, corresponding approxi- 
mately to maximum in drag with a minimum in lift 
and a minimum in drag with zero lift. Because of the 
contributions to the drag from both upper and lower 


vorticies forming behind the cylinder, the drag varies 
with double frequency compared to that of the lift. 


3.3. Wake Vortex Propagation 

The objective of this case is to investigate how 
a wingtip vortex is preserved in the present compu- 
tational procedure. Since wingtip vorticies can ex- 
ist for hundreds of chord-lengths behind the wing, a 
tip vortex generated by a large airplane can influence 
following aircraft. Blade/ vortex interaction on rotor- 
craft and tip vortex cavitation on ship propeller blades 
are also areas where an accurate tip vortex simula- 
tion has a significant role. The computational study 
of the near-field behavior of a wingtip vortex using the 
artificial compressibility method, 27,28 done in conjuc- 
tion with the experimental study, 29 indicated that it 
is possible to predict the mean flow with some degree 
of accuracy. It also has been shown that the success in 
capturing the flow features of the tip vortex depends 
on the grid resolution and the turbulence modeling. 
The static pressure coefficient in the vortex core was 
not predicted well, however, unless grid refinement in 
the core region was performed. One of the challeng- 
ing problems in vortical flow simulation is to deter- 
mine the source of numerical inaccuracy, whether it 
is from numerics, grid, or turbulence model. In or- 
der to isolate undesired grid effect from the turbulence 
modeling questions, Dacles-Mariani 27 studied wingtip 
vortex propagation using the artificial compressibility 
method. This wake vortex problem was selected as 
a validation case for the current fractional-step algo- 
rithm because of the challenging nature of the prob- 
lem. 

The schematic of the experimental test section 
and the computational domain with an H-H grid topol- 
ogy is shown in Figure 8. The wing has an aspect ratio 
of 0.75 and 10° angle of attack. The Reynolds number 
for this flow is 4.6 million based on the chord length. 
The computational domain includes the region from 
the trailing edge of the wing (x/c=1.0) to the 0.673 
of the chordlength c, downstream of the wing. Ex- 
tensive experimental data 29 is available at x/c = 1.0, 
x/c = 1.12, x/c = 1.24, x/c = 1.447, and x/c = 1.673. 
The experimental velocity profile at x/c = 1.0 sta- 
tion is used as inflow boundary conditions, and exit 
velocity components are extrapolated from the inte- 
rior. Presure distributions at boundaries are calcu- 
lated from the compatibility condition. The computa- 
tions are carried out for two grid levels; a coarse grid 
which contains 36 x 42 x 42 grid points, and a relatively 
fine grid with dimensions of 36 x 82 x 82. 
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Figure 8 : Computational domain and grid for 
wake vortex flow calculations. 

A typical convergence history for the coarse grid 
is plotted in Figure 9. The time step is set to 0.01 
for this calculation. Dashed line shows the maximum 
changes in the dependent variables U\ U v , and U' 
momentum equations from time level nton+1 (equa- 
tions 4 and 5). The solid line represents the maximum 
value of divergence of the velocity error. The solution 
is considered converged in 1000 iterations. The CPU 
time required for this computation is less than one 
hour for coarse grid and about 3.5 CRAY-C90 hours 
for finer grid. 



Figure 9: Convergence history for wake vortex 
calculation 


Initially, the computations were carried out for 
the coarse grid. It was found that the vortex core ve- 
locity peak values were underpredicted. The grid res- 
olution was inceased by doubling the number of grid 
points in k, and / direction, from grid dimension of 
36 x 42 x 42 to 36 x 82 x 82. The prediction of the peak 
values at the vortex core was improved, but not sub- 
stantially. The numerical results indicated that there 
is an excessive amount of numerical dissipation at the 
vortex core as it progresses downstream. This was 
consistent with the findings from the tip vortex study 
using artificial compressibility 27 . The excessive nu- 
merical diffusion at the vortex core was reduced when 
the production term in equation 32 is used instead of 
the magnitude of the vorticity (equation 31), as sug- 
gested by reference 27. It was observed that a numer- 
ical experimentation was necessary to find the right 
value of xf actor in equation 32. In order to remove 
this dependency, the production term in equation 33, 
which is the minimum value of the magnitude of vor- 
ticity and the norm of the strain tensor, was used. A 
very good agreement obtained between the experimen- 
tal data and computed results obtained by the last two 
approaches. Figures 10 through 14 compare computed 
results obtained from this procedure with the experi- 
mental data. 


■ 

Experiment, Zilliac & Chow (1991) 

Computed 

results 

Convective 
terms dlff.. 

Grid size 
jxkxl 

Baldwin- Bart model 
production term 

— 

3-rd upwind 

36 x 82 x 82 

nrinf |s | , |w|) 


5-th .jnu/mH 

36 x 42 x 4? 

TTmn( Is 1 . Iwl) 


5-th upwind 

36 x 82 x 82 

| w |+2.O*tran(0A| s |-| w |) 

.... 

5-th upwind 

36 x 82 x 82 

|w| 

— 

5-th upwind 

36 x 82 x 82 

mint |s | , |w|) 


Table 2 : Legend of the computed results pre- 
sented in Figures 10 through 14 


The legend for figures 10 through 14 is given in Table 
2. Symbols represent the experimental data, 29 and 
lines represent the computed results. Chain-dotted 
lines show the results from b th order differencing from 
the coarse grid with dimensions of 36 x 42 x 42. The 
production term in equation 33 was used for this case. 
Chain-dashed lines represent the results from the rela- 
tively fine grid with dimensions of 36 x 82 x 82. The or- 
der of differencing and the production term are same as 
the previous case. Solid lines show the results from 3 rd 
order upwind differencing, keeping the grid size and 
the production term unchanged. Dotted lines show 
the results from the 36 x 82 x 82 grid by using 5 th 
order differencing and using the production term in 
equation 32. The dashed lines represent the results by 
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using the production term in equation 31, which was 
originally used by Baldwin and Barth. 




Figure 10 : Axial progression of flow quantites 
along vortex coreline. 

Figure 10 shows the axial progression of flow quan- 
tities along vortex coreline. Velocity magnitude and 
static pressure coefficient Cp at the vortex core are 
plotted. The effect of the grid resolution (chain dot- 
ted lines) can be seen from these plots. Even though 
higher order of differencing and an improved produc- 
tion term in the turbulence model are used, the vortex 
core has not been resolved accurately with the coarse 
grid. When the grid resolution is increased the agree- 
ment between the numerical results and the experi- 
mental data is very good (chain-dashed lines), reduc- 
ing the error in the viscous core region to less than 2%. 
The production terms in equation 32 and 33 give al- 
most identical results (dotted and chain-dashed lines). 
However, the dashed lines clearly show the effects of 
the turbulence model. Using the magnitude of the vor- 
ticity in the production term gives excessive diffusion 
at the vortex core as it progresses downstream. The 
effect of order of differencing can be seen by comparing 
the solid lines with chain-dashed lines. The numerical 


dissipation is automatically computed in upwind dif- 
ferencing, and the amount of dissipation is large when 
third-order flux difference splitting is used at this grid 
level. In order to reduce this numerical dissipation, 
one can use finer grid or can increase the stencil in 
the upwind-biased differencing. Since the cost of in- 
creasing the accuracy of the differencing is much less 
than of increasing the grid size, the 5 th order upwind 
differencing is used for the rest of the cases. It should 
be pointed out that the overall spatial accuracy of the 
method is second-order even though 5 th order upwind 
differencing is used for convective terms. The reason 
for that is because the volume and surface area vector 
are evaluated as second-order, a simple averaging is 
used for the metric terms at half points location, and 
a second order centra] differencing is used for viscous 
terms. However, increasing the stencil in upwind dif- 
ferencing has a significant effect in reducng the amount 
of numerical dissipation, compared to lower-order dif- 
ferencing. 




Figure 11: Axial progression of the vortex core- 
line 


The vortex core location is plotted in Figure 11. 
In the coarse grid solution, the vortex core location 
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(z/c) remains unchanged because of the grid resolu- 
tion. The use of various production terms in the tur- 
bulence model has a significant effect in the vortex 
core values (Figure 10). Figure 11 shows that vortex 
core location does not change with using different pro- 
duction terms. However, grid resolution and order of 
accuracy in the differencing have influence in location 
of the vortex coreline. 


Figure 12 shows the comparison of velocity mag- 
nitude across the wake vortex at three interior sta- 
tions ( x/c = 1.12, x/c = 1.24, and x/c = 1.447) and 
at the exit boundary ( x/c = 1.673). Since the ex- 
perimental velocity profile is prescribed at the inflow 
boundary ( x/c = 1.0), the results at x/c = 1.0 are 
not plotted. Resolving the vortex core peak values 
was the focus point for these calculations. As one can 
see from Figure 8, the grid resolution for the the wind 
tunnel wall boundary layer is not sufficient, especially 
for the Reynolds number of 4.6 Million. It should also 
be noted that the turbulence Reynolds number Rt in 
the Baldwin-Barth turbulence model is set to 1 at the 
inflow boundary. In actuality, the inflow boundary for 
the computational domain lies in the wake region for 
the experiment. Therefore, the boundary layer thick- 
ness on the wind tunnel wall is larger in the computed 
results than in the experimental measurments. This 
difference results in discrepancies between computed 
results and experimental data near the wind tunnel 
wall. The differences Eire the largest at the exit plane. 
Comparison of crossflow velocity across the wake vor- 
tex at four different stations is given in Figure 13. Us- 
ing different production term in the turbulence model 
does not have a great influence in crossflow velocities 
as it has in the velocity magnitude (compare Figure 12 
and 13). Since the crossflow velocity is zero at the vor- 
tex core, the dissipation introduced by the turbulence 
model, the grid spacing and the order of accuracy di- 
rectly effects the axial velocity components. The com- 
parison of Cp acros wake vortex at the inflow bound- 
ary, at one interior station ( x/c — 1.24), and at the 
exit is plotted in Figure 14. Since the pressure has 
not been prescribed from experimental data, the com- 
puted Cp values are compared with the data at the 
inflow and the outflow boundaries. Obtaining a good 
comparison for the pressure at boundaries is especially 
encouraging for time-accurate comptations. The leg- 
end for computed results in these figures are given in 
Table 2. Dotted and chain-dashed lines show very 
good comparison with the experimental data. Find- 
ings from the previous figures about the effects of the 
grid resolution, order of accuraccy in differencing, and 
turbulence model are also valid for Cp plots. 






Figure 12: Comparison of velocity magnitude 
across wake vortex. 
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cross-flow velocity cross-flow velocity cross-flow velocity cross-flow velocity 








Figure 13: Comparison of crossflow velocity 

across wake vortex. 



y/x distance ( @ x/c = 1.673 ) 

Figure 14: Comparison of Cp across wake vor- 
tex at inflow, exit and one interior station. 


4. Conclusion 

A fractional-step algorithm for computing both 
steady-state and time-accurate solutions to the incom- 
pressible Navier-Stokes equations has been presented. 
Using the mass fluxes as dependent variables in mo- 
mentum equations automatically satisfies the incom- 
pressibility condition. The staggered grid orientation 
simplifies the pressure boundary condition in devising 
a Poisson solver and provides more compact stencils. 
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The use of upwind differencing with a non- factored im- 
plicit line-relaxation scheme in the solution of momen- 
tum equations provides fast convergence with relaxed 
stability restrictions. 

The computed results showed good comparison 
with analytical solution of a laminar Couette flow and 
with experimental and other numerical results for an 
unsteady laminar flow over a circular cylinder. The 
validated solution procedure is applied to a three- 
dimensional, steady, turbulent wingtip vortex wake 
flow. The algorithm does a very good job in resolv- 
ing the vortex core when 5 <A -order upwind differenc- 
ing and proper production term in the Baldwin-Barth 
one-equation turbulence model are used with adequate 
grid resolution. 

Since the incompressibility condition is satisfied 
automatically at every physical time-step, the current 
fractional-step solver (INS3D-FS) will be applied to 
unsteady, three-dimensional flow problems. 
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